Optimal estimates of the diffusion coefficient of a single Brownian trajectory 
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Modern developments in microscopy and image processing are revolutionizing areas of physics, 
chemistry and biology as nanoscale objects can be tracked with unprecedented accuracy. The goal 
of single particle tracking is to determine the interaction between the particle and its environment. 
The price paid for having a direct visualization of a single particle is a consequent lack of statistics. 
Here we address the optimal way of extracting diffusion constants from single trajectories for pure 
Brownian motion. It is shown that the maximum likelihood estimator is much more efficient than 
the commonly used least squares estimate. Furthermore we investigate the effect of disorder on the 
distribution of estimated diffusion constants and show that it increases the probability of observing 
estimates much smaller than the true (average) value. 
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I. INTRODUCTION 

Single particle tracking dates back to the classic study 
of Perrin on Brownian motion (BM) [l[ . It generates the 
position time series of an individual particle trajectory 
B(t) in a medium (see, e.g., Refs. 0, 3) ^-^d when prop- 
erly interpreted, the information drawn from a single, or 
a finite number of trajectories, can provide insight into 
the mechanisms and forces that drive or constrain the 
motion of the particle. The method is thus potentially a 
powerful tool to probe physical and biological processes 
at the level of a single molecule [J|. At the present time, 
single particle tracking is widely used to characterize the 
microscopic rheological properties of complex media |5[ , 
and to probe the active motion of biomolecular motors 
Q. In biological cells and complex fluids, single particle 
trajectory (SPT) methods have, in particular, become in- 
strumental in demonstrating deviations from normal BM 
of passively moving particles (see, e.g., Refs. |7l-[lG|). 

The reliability of the information drawn from SPT 
analysis, obtained at high temporal and spatial resolu- 
tion but at expense of statistical sample size is not always 
clear. Time averaged quantities associated with a given 
trajectory may be subject to large fluctuations among 
trajectories. For a wide class of anomalous diffusions de- 
scribed by continuous-time random walks, time-averages 
of certain particle's observables are, by their very nature, 
themselves random variables distinct from their ensem- 
ble averages [il| . An example is the square displacement 
time-averaged along a given trajectory, which differs from 
the ensemble averaged mean squared displacement |12| . 



By analyzing time-averaged displacements of a particular 
trajectory realization, subdiffusive motion can actually 
look normal, although with strongly diffe ring diffusion 
coefficients from one trajectory to another [l3|. 

Standard BM is a much simpler and exceedingly well- 
studied random process [ij] than anomalous diffusion, 
but still it is far of being as straightforward as one might 
be tempted to think. Even in bounded systems, despite 
the fact that the first passage time distribution has all 
moments, first passages to a given target of two inde- 
pendent identical BMs, starting at the same point in 
space, will most likely occur at two distinctly different 
time moments J15j . revealing a substantial manifesta- 
tion of sample-to-sample fluctuations. Ergodicity, that 
is, equivalence of time- and ensemble-averages of square 
displacement holds only in the infinite sample size limit. 
In practice, it means that standard fitting procedures 
applied to finite (albeit very long) trajectories of a given 
particle will unavoidably lead to fluctuating estimates of 
the diffusion coefficient D. Indeed, variations by orders 
of magnitude have been observed in SPT measurements 
of the diffusion coefficient of the Lad repressor protein 
along elongated DNA ^ (see also Section |VE|. Sig- 
nificant sample-to-sample fluctuations resulting in broad 
histograms for the value of the diffusion coefficient have 
been observed experimentally for two-dimensional (2D) 
diffusion in the plasma membrane [3| , as well as for diffu- 
sion of a single protein in the cytoplasm and nucleoplasm 



of mammalian cells [17 1 . 
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Such a broad dispersion of the value of the diffusion 
coefficient extracted from SPT measurements, raises im- 
portant questions about the correct or optimal method- 
ology that should be used to estimate D. Indeed, these 
measurements arc performed in rather complex environ- 
ments and each SPT has its own history of encounters 
with other species, defects, impurities, etc., which in- 
evitably results in rather broad histograms for observed 



D. On the other hand, it is highly desirable to have a 
reliable estimator of the diffusion coefhcient even for the 
hypothetical "pure" cases, such as, e.g., unconstrained 
standard BM. A reliable estimator should produce a dis- 
tribution of D as narrow as possible and with the most 
probable value as close as possible to the ensemble av- 
erage one. A knowledge of the distribution of such an 
estimator could provide a useful gauge to identify effects 
of the medium complexity as opposed to variations in 
the underlying thermal noise driving microscopic diffu- 
sion. Commonly used methods of extraction of D from 
the SPT data are based on a least square (LS) estimate 
of the time-averaged square displacement and some of its 
derivatives (see, e. g., [ 31, IitI. Il8j and the next section). A 
recent study, Ref. [13, focussed on estimators for D for 
ID BM, the statistics of which is amenable to analyti- 
cal analysis. Several methods for estimating D from the 
SPT data were studied and it was shown that a com- 
pletely different approach - consisting of maximizing the 
unconditional probability of observing the whole trajec- 
tory - is superior to those based on the LS minimization. 
As a matter of fact, at least in ID systems the distribu- 
tion of the maximum likelihood (ML) estimator of the 
diffusion coefficient not only appears narrower than the 
LS ones, resulting in a smaller dispersion, but also the 
most probable value of the diffusion coefficient appears 
closer to the ensemble average D [l^ . 

In this paper we focus first on the case of pure standard 
BM and calculate exactly, for arbitrary spatial dimension 
d, the distribution P{u) of the maximum likelihood esti- 
mator 



I 
T 



ta+T 



dt 



to 



B^(t) 
E{B2(i)}' 



(1) 



of the diffusion coefficient of a single BM trajectory B(t). 
The parameter ig here is the lag time (at which the mea- 
surement is started) which can be set equal to zero for 
standard BM without any lack of generality. However 
for anomalous diffusion, or BM in presence of disorder ig 
will play a significant role. The symbol E{. . .} denotes 
ensemble average, so that 



E{B2(t)} ^2dDt, 



(2) 



D being the ensemble-average diffusion coefficient. Con- 
sequently, the random variable u is defined as the ratio of 
the realization-dependent diffusion coefficient, calculated 
as the weighted time-average along a single trajectory, 
and the ensemble average diffusion coefficient. Clearly, 
E{u} = 1. 

Further on, we analyze here a useful measure of 
sample-to-sample fluctuations - the distribution function 
P{uj) of the random variable 



Ul 



Ui +U2 



(3) 



where ui and U2 are two identical independent random 
variables with the distribution P{u). Hence, the distri- 
bution P{lo) probes the likelihood of the event that the 



diffusion coefficients drawn from two different trajecto- 
ries are equal to each other. 



Finally, we discuss the effect of disorder on the dis- 
tributions P{u) and P(w) for ID BM in random me- 
dia. We consider two different models of diffusion in ID 
random environments - diffusion in presence of a ran- 
dom quenched potential with a flnite correlation length, 
as exemplified here by the Slutsky-Kardar-Mirny model 
[20|, and diffusion in a random forcing landscape - the 
so-called Sinai model l21[. The former is appropriate 
for diffusion of proteins on DNA, which is affected by the 
base-pair reading interaction and thus is sequence depen- 
dent, while the latter describes, for example, the dynam- 
ics of the helix-coil boundary in a melting heteropolymer 
[23 |. Note that in the former case, at sufficiently large 
times, one observes a diffusive- type motion with B^ ^ f, 
while in the latter case dynamics is strongly anomalous 
so that Bt is logarithmically confined, B^ '^In t. 



The paper is outlined as follows: In Section [IT] we re- 
call some common fitting procedures used to calculate the 
diffusion coefficient from single particle tracking data. In 
Section [nil we focus on the maximum likelihood estima- 
tor and, generalizing the approach developed in Ref. |19j 
for ID systems, obtain new results for the moment gener- 
ating function $(cr) and the probability density function 
P{u) of the ML estimator for arbitrary spatial dimension 
d. In that Section we also obtain the asymptotical be- 
havior of the probability distribution function P{u), as 
well as its kurtosis and skewness. Next, in Section IIVI 
we focus on the probability distribution function of the 
random variable w - a novel statistical diagnostics of the 
broadness of the parental distribution P{u) which probes 
the likelihood of the event that two estimates of the diffu- 
sion coefficient drawn from two different trajectories are 
the same. Further on, Section|V]pi'esents a comparison of 
the commonly used least squares estimator and the max- 
imum likelihood estimator. We show that the latter out- 
performs the former in any spatial dimension d producing 
a lower variance and the most probable value being closer 
to the ensemble average value. Next, in Section I VII we 
focus on Brownian motion in presence of disorder. As 
exemplified by two models of dynamics in systems with 
quenched disorder - Sinai diffusion (random force) and 
Slutsky-Kardar-Mirny model (random potential), disor- 
der substantially enhances the importance of sample-to- 
samplc fluctuations. We show that the observation of 
values of the diffusion coefficient significantly lower than 
the ensemble average becomes more probable. We show, 
as well, that as the strength of disorder is increased, the 
distribution P{lj) undergoes a surprising shape-reversal 
transition from a bell-shaped unimodal to a bimodal form 
with a local minimum at w = 1/2. Finally, we conclude 
in Section IVIII with a brief recapitulation of our results 
and some outline of our further research. 



II. FITS FOR THE DIFFUSION COEFFICIENT 
OF A SINGLE TRAJECTORY 

To set up the scene, we first briefly recall several fitting 
procedures commonly used to calculate the diffusion co- 
efficient from the SPT data. More detailed discussion can 
be found in Refs. |a. ll7Hl9| . We focus here on estimators 
which yield a first power of D. Non-linear estimators, e.g. 
a mean maximal excursion method |23| which has been 
used to study anomalous diffusion and produces vZ3, will 
be analyzed elsewhere. 

One of the simplest methods consists in calculating a 
least squares estimate based on the minimization of the 
integral 



where the trajectory B(t) is appropriately discrctizcd. 
Differentiating the logarithm of L with respect to D and 
setting d\n.L/dD = 0, one finds the maximum likelihood 
estimate of D, which upon a proper normalization is de- 
fined by Eq. ([T]). 

Below we will derive the distribution function P{u) of 
the ML estimator u and compare it against numerical 
results for the distribution function of the LS estimator 
Us ior d = \ , 2,3. 



III. DISTRIBUTION OF THE ML ESTIMATOR 



A. The moment generating function 



dt{B^{t)-i{t)y 



(4) 



where the diffusion law l{t) is taken either as a linear, 
l(t) = 2dDit, or an affine function, l{t) = 2dDat + ha- In 
particular, for the linear case the least squares minimiza- 
tion yields the following linear-least-squares estimator: 
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(5) 



where A is the normalization factor, A = 3/2dDT^, con- 
veniently chosen so that E{u;s} = 1. 

A second, more sophisticated, approach is based on 
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(6) 



which is the temporal moving average over a sufficiently 
long trajectory B(t) produced by the underlying process 
of duration T ^ t. The diffusion coefficient is then ex- 
tracted from fits of 6'^{t), or from a related least squares 
estimator, which is given by the following functional of 
the trajectory 



us = 7^ 



dttSr 



it), 



(7) 



where A is the same normalization constant as in Eq. ([5|) . 
Note that the random variable us is again conveniently 
normalized so that Eju^} = 1, which enables a direct 
comparison of the respective distributions of different 
estimators. As shown in Ref. [19|, us only provides a 
slightly better estimate of D than uis- A conceptually 
different fitting procedure has been discussed in Ref. [I3 
which amounts to maximizing the unconditional proba- 
bility of observing the whole trajectory B(t), assuming 
that it is drawn from a Brownian process with mean- 
square displacement 2dDt (see Eq. ^). This is the max- 
imum likelihood estimate which takes the value of D that 
maximizes the likelihood of B(t), defined as: 



L = '[l{A'n:dDt)~'^^^exp 



t=o 



AdDt 



(8) 



Let $(cr) denote the moment generating function of 
the random variable u defined in Eq. ([1]) , 



$(cr) =E{e-'""}. 



(9) 



The squared distance from the origin, of d-dimcnsional 
BM at time t for a given realization, decomposes into the 
sum 



B\t)^Y.Bfit), 



(10) 



i=l 



Bi{t) being realizations of trajectories of independent ID 
BMs (for each spatial direction). Thus, $(0-) factorizcs 



$(a) = G'^(a), 



where 



G{a) = E <^ cxp 
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Here, in order to calculate G{a), we follow the strategy 
of Ref. [101 and introduce an auxiliary functional: 



*(a;,i) =E^ <^ cxp 
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(13) 



where the expectation is for a BM starting at x at time 
t. We derive a Feynman-Kac type formula for ^(a;,t) 
considering how the functional in Eq. (|T3| evolves in the 
time interval (i,i -I- dt). During this interval the BM 
moves from x to x + dBi{t), where dBi[t) is an infinitesi- 
mal Brownian increment such that EdB{(ii?i(t)} = and 
^dB{dBf{t)} = 2Ddt, where E^b denotes now averag- 
ing with respect to the increment dBiit). For such an 
evolution we have to order dt: 
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Expanding the right-hand-side of the latter equation to 
second order in dB^ (t) , hnear order in dt and performing 
averaging, we find eventually the following Schrodinger 
equation: 



di 



= -D 



d'^^{x,t) GX^ 



dx^ 



2dDTt 



^{x,t). (15) 



The solution of this equation has been obtained in 
Ref. [13 and gives 

G(a)=vi;(o,0)=(/o(v/8^))"'^', (16) 



where /o(.) is the modified Bessel function [2J1. Conse- 
quently, we find the following general result 



$(a)=(/o(v/8^)) 



-d/2 



(17) 



Note that $((t) is independent of T and D, as it should 
in virtue of the scaling properties of the BM. 
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FIG. 1: (Color online) The distribution P{u) in Eq. ^ for 
ML estimates of BM: the blue solid line (left) corresponds to 
d = 1, the red line (middle) to d = 2 and the green line (right) 
to rf = 3. The open circles present the results of numerical 
simulations. In the inset the dashed lines correspond to the 
small- ?i and large-u asymptotics in Eqs. (|26|) and (|28|) . 



B. The distribution function 

We turn next to the analysis of the distribution of the 
ML estimator defined in Eq. ([T]). First of all, we calculate 
several first moments of u by merely differentiating the 
resuh in Eq. P71) : 
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(18) 



Consequently, one may expect that all moments tend to 
1 as d — J- 00, so that P{u) -^ S{u — 1). For fixed d, 
the variance E{m^} — E{u}" = 1/d, the coefficient of 
asymmetry 7a = 8/3\/d and the kurtosis 7e = ll/rf. All 
these characteristics vanish when d — > 00. 

Note next that since Io{y/8a-/d) = Jo{i\/8(7/d), the 
poles of ^{(j) are located at cr = -~d^l/8, where 7^ is the 
fcth zero of the Bessel function Jo(.) [24| . Consequently, 
for even d, we can straightforwardly find P{u) in form 
of an infinite series in the zeros of the Bessel function 
Jo(.). For d = 2, ^{a) has only simple poles so that the 
expansion theorem gives 
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For d = 4 and d ~ 6 the standard residue calculus yields 
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where 

ak = 97|jf (7fc)w'-367,'J?(7fc)ii+4J|(7fe)(7fe-8). (22) 

Similar results can be readily obtained for greater even 
d. 

For arbitrary d, including odd values, the distribution 
P{u) is defined by inverting the Laplace transform and 
is given by the following integral: 
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(24) 
ber(a;) and bei(a;) being the 0th order Kelvin functions 
11. 

Finally, we consider the small-u and large-u asymp- 
totic behavior of the probability density function P{u). 
To extract the small-u asymptotic behavior of P{u) we 
consider the large-u form of $(cr). From Eq. ([T7|) we get 
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(20) &s a ^f 00. Consequently, we find the following strongly 



non-analytical behavior: 

D/ ^ tA ^d/4 rd ( d\ I d + 2 

(26) 
The large-u behavior of the distribution P{u) is defined 
by the behavior of the moment generating function $(cr) 
in the vicinity of cr* = —djf/S, 
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Consequently, we find that for m — > cx3. P{u) decays as 
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... (28) 
This behavior is, of course, consistent with the series ex- 
pansions in Eqs. (|T9|) . (|20|) and ((2T|) . Our results on the 
distribution P{u) are summarized in Fig. [TJ 




value? Of course the distributions and thus moments 
of these two random variables are the same, however a 
measure of their relative dispersion can be deduced by 
studying the distribution function P(w) of the random 
variable w [H, ll^, defined in Eq. ([S])- This distribution 
is given explicitly by [26| 
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(29) 



and hence, it suffices to know P{u) in order to determine 
P(c^). 

For d = 2, (and, in fact, for any other even d), P{uj) 
can be evaluated exactly. Plugging Eq. (|19p into ([2^ . we 
get: 
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(30) 



Performing the sum over /, we arrive at the following 
result for the distribution P{uj) in 2D systems 
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Numerically obtained distributions P(w) for d = 1, 2 and 
d = 3 dimensional systems are presented in Fig. [2] No- 
tice that in all dimensions u ~ 1/2 is the most proba- 
ble value of the distribution P{uj) so that most proba- 
bly ui = U2. Nevertheless, the distributions are rather 
broad which signifies that sample-to-samplc fluctuations 
are rather important. 



COMPARISON OF THE LS AND ML 
ESTIMATORS. 



FIG. 2: (Color online) Distribution P{cu) for standard BM 
obtained from Eqs. (|23|) and (|29|) . P{ijj = 1/2) increases 
with the dimension d. The blue (lower) solid line corresponds 
to d = 1, the red (middle) line to d = 2 and the green (upper) 
line to d = 3. Open circles with the same color code (and 
relative positions) present the results of numerical simulations 
for the LS estimator us defined in Eq. ((T]). Note that an 
apparent coincidence of the results for the distributions P(tj) 
for the ML estimator in 2D and that for the LS estimator us in 
3D is accidental. It just signifies that the former outperforms 
the latter. 



IV. THE DISTRIBUTION OF THE RANDOM 
VARIABLE Lj 

Suppose next that we have two different independent 
realizations of BM trajectories, Bi(i) and B2(t) which we 
use to generate to independent random variables ui and 
U2. A natural question arising about their suitability as 
estimators is how likely is it that they will have the same 



We will now show that the ML estimator defined in 
Eq. ([T]) substantially outperforms the MS estimator as 
defined in Eq. ([7]) in any spatial dimension d. This is a 
very surprising result as one would intuitively expect, and 
it is often stated in the literature, that using the process 
6t has the effect of a reducing the fluctuations of the 
estimate of D because the process is partially averaged 
in time. 

To demonstrate this, we present in Fig.[3]a comparison 
of the analytical results for P{u) of the ML estimator 
with the corresponding distributions of the LS estimator 
for Us obtained numerically. 

Indeed, we find that the variance of the distribution 
P{us) equals 1.38, 0.66 and 0.44 for d = 1,2 and 3, re- 
spectively. The distribution of the ML estimator appears 
to be substantially narrower so that the variance is sig- 
nificantly lower, 1, 1/2 and 1/3. Moreover, the most 
probable values of u are closer to the ensemble average 
value E{u} = 1 than the most probable values of P{us) 
to Ejwa} = 1: we observe that the distribution P{us) at- 
tains its maximal values at us ~ 0.15, 0.33 and 0.47 for 
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FIG. 3: (Color online) Comparison of P{u) in Eq. ((23l) (solid 
lines) and the results of the numerical simulations for the 
distribution P{us) (histograms). From left to right: d = 1, 2 
and 3. 



d = 1,2 and 3, respectively, while the corresponding max- 
ima of the distribution P{u) are located at m « 0.28, 0.47 
and 0.6. Last but not least, the distribution P(ujs) ap- 
pears to be significantly broader than P(w), as revealed 
by Fig. [H The worst performance of the LS estimator us 
is in ID systems in which the distribution P{ios) has a bi- 
modal Af -shape with a local minimum at cog = 1/2, and 
maxima (most likely values) around 0.1 and 0.9. This 
means that the values ui and U2 drawn from two dif- 
ferent trajectories will most probably be different by an 
order of magnitude! 



VI. ID BROWNIAN MOTION IN PRESENCE 
OF DISORDER. 

In this final section we address the question of how 
the distribution P{u) of the ML estimator of a single 
trajectory diffusion coefficient will change in presence of 
quenched disorder. We will consider two different models 
of BM in random ID environments - diffusion in presence 
of a random correlated potential and diffusion in presence 
of a random force. 



A. Diffusion in presence of a random potential 

First we consider a BM in a ID inhomogcncous en- 
ergy landscape, where disorder is correlated over a fi- 
nite length ^c- This model gives a simple description 
of diffusion of a protein along a DNA sequence, for in- 
stance, where the particle interacts with several neigh- 
boring base pairs at a time [20] ■ The total binding en- 
ergy of the protein is assumed to be a random variable. 
When the particle hops one neighboring base further to 
the right or to the left, its new energy is highly cor- 
related to the value it had before the jump. Slutsky 
et al. [20| modeled this process as a point-like parti- 
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FIG. 4: (Color online) Distribution P{u) for a particle diffus- 
ing in a random energy landscape with variance a^, correla- 
tion length ^c = 20 and various disorder strengths e — /3a. 
From right to left: blue histogram corresponds to e = 0.5, red 

- to e = 0.8, green - to e = 1, yellow - to e = 2 and brown 

- to e = 3. The walk duration is T = 10^, io = and av- 
erages are taken over 12, 000 walks occurring in independent 
landscapes. For comparison we present the distribution (solid 
black curve) for standard ID BM {e = 0). The correspond- 
ing distributions P(tj) are shown in the inset {P{uj — 1/2) 
decreases with increasing e). 



clc diffusing on a ID lattice of unit spacing with ran- 
dom site energies {Ui}, whose distribution is Gaussian 
with zero mean, variance a'^ and is correlated in space as 
{{Ui - Ujf) = 2cr2[l - exp(-|i - Jl/'Cc)]. At each time 
step, the particle located at some site i jumps to the left 
or to the right with probabilities pi a exp[/3(/7i — C/^i-i)] 
and qi ex exp[/3(t/i — C/i+i)], respectively, where Pi + Qi = 
1. Diffusion is asymptotically normal for any disorder 
strength e = /3a. Nevertheless, the particle can be 
trapped in local energy minima for long periods of time. 
During an extended intermediate time regime, it is ob- 
served that first passage properties fluctuate widely from 
one sample to another (20|. 

Our numerical simulations reveal that disorder has a 
dramatic effect on the distributions P{u) and P{oj)- As 
shown by Fig. |4l the distribution P{u) broadens signif- 
icantly in the small u regime: very small values of the 
time-average diffusion constant (compared to the thermal 
and disorder average) become increasingly more probable 
as the disorder strength increases. However, the right tail 
of P{u) is much less affected. Similarly, two independent 
measurements arc likely to differ significantly, even in 
moderately disordered media (sec inset of Fig. 01). When 
e ~ 0.8, the distribution P{uj) undergoes a continuous 
shape reversal transition - from a unimodal bell-shaped 
form to a bimodal A/-shape one with the minimum at 
cj = 1/2 and two maxima approaching and 1 at larger 
disorder strengths. Unfortunately, it does not seem possi- 
ble to obtain this critical value analytically. Even for the 
case of a pure Brownian motion considered in Ref. |la | , 
such an analysis appears to be extremely difficult. 



Therefore, for e > 0.8 samplc-to-samplc fluctuatfons 
becomes essential and it is most likely that the diffusion 
coefficients drawn from two different trajectories will be 
different. 



P{u) 




In u 
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FIG. 5: (Color online) In panel (a), distribution P{u) for 
Sinai diffusion and different strengths of the disorder e. From 
right to left: blue histogram corresponds to e = 0.05, red - to 
e = O.I, green - to e = 0.2 and yellow - to e = 0.3. The solid 
black curve depicts the distribution P(u) for ID BM (e = 0), 
and the corresponding distributions P{uj) are shown in the 
inset with the same color code (for ui close to or 1, P{u]) 
increases with e). In panel (b), distribution P{u]) for Sinai 
diffusion with e = 0.1, the integration time is set to T = 10^ 
and to is varied. P{u} = 1/2) decreases with increasing to: the 
different curves corresponds to (from darker gray to lighter 
gray): to = 5 (blue), to = 5 x 10^ (red), to = 5 x 10^ (green) 
and to = 5 X 10^ (yellow). In panel (c), the lag time is set 
to to = 5 and the integration time T is varied. P{u} = 1/2) 
decreases with increasing T: 10^ (blue), 5x10'' (red), 1, 5x 10'' 
(green) and 2, 5 x 10^ (yellow). 



variables with distribution P{ei) = \5{ei — e) + \5{ei + e) 
and the strength of disorder e is bounded away from and 
1 . It is well-known that in the large-t limit the model pro- 
duces an anomalously slow sub-diffusion (x^(<)) ^ In (t), 
where the angle brackets denote averaging with respect 
to different realizations of disorder. At shorter times, 
however, one observes an extended stage with a tran- 
sient behavior which is substantially different from the 
asymptotic one. As a consequence, the statistics of w, 
defined in equation ([T]), and w will depend not only on 
the integration time T but also on the lag time to from 
which a single trajectory is analyzed. 

Wc have numerically computed the distributions P(u) 
and P(w). In Fig. [Sji we present the dependence of the 
u and a; statistics on the strength of the disorder e. As 
in the previous disordered potential case, we find that 
the maximum of P{u) shifts toward zero as the disor- 
der gets stronger. For comparison, the solid black line in 
Fig-Eb represents P{u) observed for standard BM, e = 0. 
Moreover, the stronger the disorder is, the broader the 
distribution P{u) becomes, yielding more peaked max- 
ima in P{ijj) (see the inset of Fig. [SJi). We also note that 
P{oj) has a bimodal M-shapcd form even for the weakest 
disorder e = 0.005 that we have considered, suggesting 
that the zero-disorder limit is non-analytic compared to 
the continuous transition observed for diffusion in the 
random energy landscape of the previous section. 

In Fig. [nt and Fig. [S]:, we show the statistics of cj for 
different values of to a-nd T . Increasing to (or T) wc ob- 
serve that P{uj) changes from an almost uniform form, 
for which any relation between ui and U2 is equally prob- 
able, to a bimodal M-shaped distribution, which signifies 
that in this regime two ML estimates ui and U2 will most 
likely have different values. Bimodality is a property of 
the Sinai regime [la], as also noticed earlier in Ref. j27| . 
and thus it shows up only at sufficiently long times when 
the trajectories follow the asymptotic ultra-slow diffu- 
sion. The distribution P{ui) is remarkably sensitive to 
the characteristic aging of Sinai diffusion. 

As a final observation, one may also study the statistics 
of u for trajectories evolving in the same random force 
field. In this case (not shown here) one gets a narrow 
distribution for of u and a unimodal P{lo) that converges 
to a delta singularity at a; = 1/2 when the disorder be- 
comes infinite. This is due to the fact, known as Golosov 
phenomenon, that two trajectories in the same disorder 
will move together |28j . 



B. Diffusion in presence of a random force 

We discuss now the effect of disorder on the distribu- 
tions P{u) and P{oj) for ID BM in presence of a quenched 
uncorrelated random force - the so-called Sinai diffusion 
(2l|. In this model, one considers a random walk on a 
ID infinite lattice and site i dependent hopping probabil- 
ities: Pi ~ \ — Ei for hopping from i to the site i + 1 and 
Qi ^ ■^ + £i for hopping to the site i — 1. Here, Ei are in- 
dependent, uncorrelated, identically distributed random 



VII. DISCUSSION 

We have analyzed the reliability of the ML estimator 
for the diffusion constant of standard Brownian motion 
and shown its superiority over the more commonly used 
LS estimator in a number of important aspects, notably 
the variance of the estimator, the proximity of the most 
probable value to the true mean value and the distri- 
bution of the random variable a; which is a measure of 



the extent to which two estimations of D vary. Going be- 
yond the important test case of pure Brownian motion we 
have also analyzed the effect of quenched disorder, model- 
ing fluctuations of the local energy landscape and forces. 
As one may have intuitively expected, the presence of 
short range disorder tends to broaden the distribution of 
the so measured value of D, as it presents an additional 
source of fluctuation. However in the Sinai model, in the 
same realization of the force field, trajectories are disor- 
der dominated and are almost independent of the ther- 
mal noise, leading to highly peaked distributions of D. 
Analytically understanding the distribution of D in the 
presence of disorder presents an interesting mathematical 
challenge which will involve analysis of the correspond- 
ing Schrodinger equation with a random drift term. Fur- 
ther interesting questions remain to be addressed, in par- 
ticular can one use the two time correlation function of 



measured trajectories to obtain better estimators ? Sin- 
gle particle tracking technology will undoubtedly further 
improve in the coming years and many interesting math- 
ematical, statistical and physical challenges will arise in 
the ultimate goal of getting the most out of the trajecto- 
ries so obtained. 
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